function [x, fval, output] = qppdbox(H, f, lb, ub, opts)
%function [x, fval, output] = qppdbox(H, f, lb, ub, opts)
%
% Quadratic programming
%   min  0.5*x'*H*x + f'*x
%   s.t. lb<=x<=ub,
% where H is symmetric positive definite.
%
% opts has the following fields:
%   maxiter: maximum number of iterations. default: 1000
%   tol: tolerance. default: 1e-6
%
% output has the following fields:
%   iter: number of iterations
%   grad_zero_part: norm of the vector grad(lb<x<ub)
%   grad_pos_part: counts how many grad(x==lb)>=0 are violated
%   grad_neg_part: counts how many grad(x==ub)<=0 are violated
%   for x to be the global min, the above three should all be zero.
%
% This program is based on the two-metric projection method.
%
% See also QUADPROG.

% Copyright 2008, Jie Chen.
% This program is free software; you can redistribute and/or
% modify it for NON-COMMERCIAL purposes. This program is
% distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY, including that of MERCHANTABILITY or FITNESS FOR A
% PARTICULAR PURPOSE.
% $Date: 2008/11/04 14:50:52$

% Bug fix:
% change ind = find(x==0&grad>0);
% to ind = find(x==lb&grad>0 | x==ub&grad<0);
% $Date: 2011/07/08 09:51:10$


if ~exist('opts','var'); opts = struct; end
if ~isfield(opts,'maxiter'); opts.maxiter = 1000; end
if ~isfield(opts,'tol'); opts.tol = 1e-6; end

n = size(f,1);
if numel(lb)==1; lb = repmat(lb,n,1); end
if numel(ub)==1; ub = repmat(ub,n,1); end

x = (ub-lb).*rand(n,1) + lb;
x_old = x;

for i = 1:opts.maxiter
  grad = H*x + f;
  
  ind = find((x==lb&grad>0) | (x==ub&grad<0));
  HH = H;
  HH(ind,:) = 0;
  HH(:,ind) = 0;
  diagind = sub2ind(size(HH), ind, ind);
  HH(diagind) = H(diagind);
  
  d = -HH\grad;
  alpha = - d'*grad / (d'*H*d);
  x = x + alpha * d;
  x = min(max(x, lb), ub);
  
  if norm(x-x_old) < opts.tol
    break;
  else
    x_old = x;
  end
end

output.iter = i;
if i == opts.maxiter
  fprintf(1,'maxiter reached.\n');
end
fval = 0.5*x'*H*x+f'*x;

%% check optimality condition
grad = H*x+f;
output.grad_zero_part = norm(grad(lb<x&x<ub));
output.grad_pos_part = length(find(grad(x==lb)<0));
output.grad_neg_part = length(find(grad(x==ub)>0));
